Continuous corrections to the molecular Kohn-Sham gap and virtual orbitals 
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We use projector operators to correct the Kohn-Sham Hamiltonian of density functional theory 
(KS-DFT) so that the resulting mean-field scheme yields, in finite systems, virtual orbitals and 
energy gaps in better agreement with those predicted by quasiparticle theory. The proposed correc- 
tion term is a scissors-like operator of the form (7 — p)5H{j — p), where I is the identity operator, 
p the density matrix of the N-particle system and 8H is either the difference between the N+l- 
and N-particle Kohn-Sham Hamiltonians or a non-self-consistent approximation to it. Such a term 
replaces the Kohn-Sham virtual orbitals of the iV-particle system by the HOMO and virtual or- 
bitals of the system with N + 1 particles in an attempt to mimic a true quasiparticle spectrum. 
Using a local density approximation (LDA) we compute the gaps of a variety of small molecules 
finding good agreement with experiment and computationally more demanding methods. For these 
systems we examine the physical origin of this gap correction and show that so-called band gap 
discontinuity, A xc , contains electrostatic contributions that do not originate from the discontinu- 
ity in the exchange-correlation potential. The similarity between the corrected and Hartree-Fock 
virtual orbitals is illustrated and the extent to which the bare LDA virtual orbitals are improved 
is considered. The lack of band-gap discontinuity and the presence of self-interaction errors in the 
proposed correction are also discussed. 
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The accurate determination of the fundamental en- 
ergy gap in many-electron systems, such as molecules 
and solids, by means of a computationally inexpensive 
mean-field method remains an open problem despite of 
recent progres a 1 ' 2 ' 3 ' ' 5 ' 6 ' , The calculation of fundamen- 
tal energy gaps based on the relevant eigenvalues of the 
Kohn-Sham Hamiltonian of Density Functional Theory 
(KS-DFT)-^ severely underestimates the experimental 
energy gap both in finite and extended system a 10 ' 11 ' 12 ' 13 . 
KS-DFT has been extremely successful in the compu- 
tation of ground state averages for which accurate den- 
sity functionals are known, such as the total energy and 
magnetic momen t 14 ' 15 ' 16 . This success is also due to the 
extremely good ratio between the accuracy and the com- 
putational demands of this method; KS-DFT not only 
yields accurate results for these quantities, but it is also 
very fast and simple to implement for most available 
exchange-correlation (XC) functionals. 

Unfortunately KS-DFT is not designed to give addition 
and removal energies in a single calculation: the Kohn- 
Sham eigenvalues have no direct physical interpretation 
except for the ones associated with the highest occupied 
molecular orbital (HOMO) 1 -^; none of the Kohn-Sham 
orbitals has a known physical meaning; the Kohn-Sham 
virtual orbitals do not enter the self-consistent cycle and 
play no role in the theory; and the Kohn-Sham Hamilto- 
nian is unphysical as we will show clearly below. 

In KS-DFT the ground state N-particle electron den- 
sity is obtained by the self-consistent solution of the 
Kohn-Sham equations, which are characterized by a non- 
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where the effective Kohn-Sham potential, v s , is a func- 
tional of the electron density of the TV-particle system, 
un, and includes the external and Hartree potentials 
as well as the XC potential, which remains unkown 
and is notoriously difficult to approximate within the 
realm of KS-DFT since it is a discontinuous functional 
of the densit y 11 ' 18 and must remove the fictitious self- 
interaction terms present in the Hartree term. 

Note that v s in Eq. |T]) is the same for all of the Kohn- 
Sham orbitals, whether occupied or not. This fact is 
in clear contrast with the Hartree-Fock (HF) equations 
which in the basis of eigenstates of the HF Hamiltonian, 
ipi, can be written as 
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where ho(r) accounts for the kinetic energy and ex- 
ternal potential, Vu(r) = J d 3 r'\r - r'l -1 ^^')! 2 and 
V i3 {v) = J dV^rOlr-rrV^r')- Since S id = 
when tpi is unoccupied the electrostatic fields that enter 
the HF equations are generated by N particles while they 
are generated by N — 1 particles when ipi is occupied 7 . 
This is precisely what one would expect from simple elec- 
trostatic considerations of a charged particle added to or 
removed from a dielectric or metallic system. While the 
HF approximation has many drawbacks associated with 
the lack of dynamical correlation, this one physically ap- 
pealling feature is totally absent in KS-DFT. Note that in 
many-body perturbation theory, the one-particle Green's 
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function in the frequency representation has poles at the 
excitation energies of the N + 1 (TV — 1) electron system 
for unoccupied (occupied) state s 19 ' 20 ! 21 . 

This is one reason, other than the discontinuity in the 
XC potential, why exact KS-DFT underestimates the ex- 
perimental gap in finite systems, while reproducing exact 
electron density as well as the total and HOMO energies. 
Consider the concrete application to finite and spheri- 
cally symmetric systems such as atoms. The exact KS- 
DFT effective potential is self-interaction free and decays 
asymptotically into the vacuum a s 5 ' 22 cx — 1/r, where r 
is the radial distace to the atomic nucleus, both for un- 
occupied and occupied states. However, in quasiparticle 
theory this potential should decay a o 23 ' 24 oc —1/r 4 for 
unoccupied states (particle added to the system). The 
oc —1/r 4 tail results from the interaction between the 
added electronic charge and the polarization charge den- 
sity it induces in the spherical symmetric system. Since 
-1/r 4 is more confining than -1/r, the KS-DFT LUMO 
energy underestimates the experimental one 3 ^, contribut- 
ing to the underestimation of the gap. These arguments 
are purely electrostatic in nature, involving only contin- 
uous functionals of the density that vanish with vanish- 
ing fractional particle numbers or in infinite systems. In 
other words the effects discussed above are proportional 
to the added charge and vanish with zero added charge 
and thus are not directly related to the well-known dis- 
continuity in the XC potential. To summarize: the ef- 
fective potential of KS-DFT makes no distinction be- 
tween occupied and virtual orbitals and, as a result, the 
screened electrostatic repulsion of the added electron by 
all the other electrons is missing in KS-DFT. 

This band-gap problem also arises in KS-DFT-based 
calculations of electronic transport through molecular 
conductor s 7 ' 37 . In Ref.[37]it is shown that exact-exchange 
DFT predicts conductances that are one order of mag- 
nitude smaller than those obtained by the Hartree-Fock 
approximation, despite of the fact that the electronic den- 
sities obtained by both methods are the same. An ad- 
ditional important problem -not addressed here- in the 
context of quantum electronic transport is the presence 
of self-interaction errors that arise through the use of ap- 
proximated XC functionals 3 ^ 3 - 7 -. In Refs. the im- 
portance of the single-particle orbitals is also highlighted. 
Atomistic simulations of electronic transport are compu- 
tationally demanding and a single-particle theory which 
is computationally inexpensive, free of self-interaction er- 
rors, able to accurately predict the molecular HOMO- 
LUMO gap and -somehow- able to provide good single- 
particle orbitals is needed. 

The rest of this paper is organized as follows: in the 
next section we discuss the correction proposed by Ce- 
hovin et ali 7 -, which should correct the above-mentioned 
deficiency of KS-DFT; in Section III we present our nu- 
merical results for the HOMO-LUMO gaps and virtual 
orbitals of a few selected molecules. We end with the 
conclusions in Section IV. 



II. CORRECTION TO THE KOHN-SHAM 
HAMILTONIAN 

These arguments show that to improve the gap val- 
ues, the iV-particle Kohn-Sham Hamiltonian needs to be 
modified in order to include information about singly- 
charged state with TV + 1 electrons. We note that the 
so-called ASCF method, where the gap is calculated 
by means of three total energy calculations as E g — 
En+i + -E^v— l — 2i?jv, is known to provide better gaps 
than those obtained directly from the Kohn-Sham eigen- 
values. Since, for all TV, the electron affinity (EA) of the 
TV-particle system and the ionization potential (IP) of 
the TV + 1-particle system satisfy: 

EA(N) = IP(N + 1) , (3) 

the fundamental gap can be obtained by means of two 
calculations of the Kohn-Sham HOMO, one for the 
N-particle system and another for the (N+l)-particle 
systemiii 2 ^ 3 -^. 

To account for these differences between the Kohn- 
Sham and quasiparticle Hamiltonians, here we apply a 
variant of the method of improved virtual orbital s 25 i 26 '' 27 . 
Given the density matrix of the TV-electron system, ,5 at, 
we define a projector over the unoccupied states of the N- 
electron system, I—pN and a mean-field-likc Hamiltonian 

H N = H* s + (I - p N )SH(I - PN ) (4) 

where SH is defined below. Because of the action of 
I ~ Pn , SH acts only on the Kohn-Sham virtual orbitals 
of the TV-particle system, leaving the ground state prop- 
erties invariant. The occupied orbitals of Eq. ^ are just 
the occupied Kohn-Sham orbitals. The virtual orbitals 
depend on how SH is defined. In applications of the so- 
called scissors operator method 2 ^, SH = el, where e is 
typically a free parameter. Following Refs. [TUIHH, 
can be defined as e = e# +1 — > anc ^ ^ s ac tion is just to 
shift the virtual-orbital energies of the TV-clcctron system 
so that Eq.© is satisfied. Here e^| i 1 ) (N) is the HOMO 
(LUMO) energy of the Kohn-Sham N+l (N) particle sys- 
tem. Such an approach corrects the virtual orbital ener- 
gies but leaves the orbitals themselves unchanged. Hence 
to obtain a mean-field-like approximation to the unoccu- 
pied quasiparticle orbitals, we define SH as 

SH = -ff^+i — H^ S . (5) 

The above equation takes into account the fact that the 
electromagnetic fields acting on accupied and unoccupied 
orbitals are generated by different numbers of particles. 
Eqs.([l][5]) replaces the virtual orbitals of the TV-particle 
system by an approximation to the HOMO and virtual 
orbitals of the TV+1 particle system, hence accounting for 
repulsion of the added particle by all the other particles, 
as well as some approximated form of screening given by 
the approximation used for XC potential in the Kohn- 
Sham Hamiltonians in Eq.©. 
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III. NUMERICAL RESULTS 

A. Computational Details 

We calculate the HOMO-LUMO gap and LUMO's of 
several isolated molecules, by means of Eqs.(|3][5]). In our 
calculations we approximate the XC potential by means 
of a spin-unpolarized LDA2£, where the correlation part 
is given as in Ref. H^. We use either the self-consistent 
(SCF) LDA density matrix of the N + 1 particle sys- 
tem or a non-self-consistent approximation (NSCF) to 
it, p N+1 w p N + |Vif }{i>L I, where V>l is the LUMO of 
the iV-particle system. We use an open source quantum 
chemistry package^ which we have tested by compar- 
ing the results of our calculations (HOMO energies and 
HOMO-LUMO gaps) with those reported in Ref. |H| as 
well as independent calculations using output from the 
GAMESS code^. From these experiences we estimate 
our numerical error bar in the orbital energies to be ~ 
1-2 mHa for LDA-based calculations and typically less 
than 1 mHa for HF calculations. We use the gaussian 
basis set 6-31G** and perform our calculations at the 
GGA geometries given in Ref. [U The geometry of Na^ 
was taken from Ref. [H. 



B. Gaps 

Calculated HOMO-LUMO gaps using different meth- 
ods are presented in Table Q] which clearly demonstrates 
that both the SCF and NSCF methods greatly improve 
the value of the gap, which is much closer to both the HF 
and ASCF gaps relative to the bare LDA gaps, in line 
with previous calculations that use essentially the same 
method±££4 In most cases both the SCF and NSCF 
methods yield gaps that are somewhat smaller than the 
HF ones. In addition the SCF is shown to produce gaps 
that are smaller than those of the NSCF. It is worth 
pointing out that the NSCF method is essentially as com- 
putationally expensive as the LDA, requiring only one 
extra iteration for the (N+l)-particle system and some 
extra post-processing, yet the gaps are greatly improved. 

In Fig[T]we consider separately the Hartree, and LDA 
exchange and correlation contributions to the SCF gap 
correction for three of these molecules. The Hartree 
contribution is obtained from the difference between the 
Hartree potentials, Vh, of the N + 1 and N particle sys- 
tem, by means of Eqs.(0][n]) with SH = Vh[f>N+i] — Vh\fiN], 
and similarly for the exchange and correlation LDA po- 
tentials. For the systems considered the main effect is 
associated to the purely electrostatic correction which in- 
creases the LDA gap by a large ammount (~ 7-8 eV for 
the molecules shown), overestimating both the Hartcc- 
Fock and experimental gaps. Both the exchange and cor- 
relation contributions to the correction given by Eqs.(T4}- 
[5]) lower the value of the gap (by about 3-4 eV for the 
molecules shown) , improve the agreement with the exper- 



Molecule 


LDA 


HF 


ASCF 


SCF 


NSCF 


H 2 


12.2 


22.5 


23.6 


20.6 


21.2 


Li 2 


1.4 


5.2 


5.6 


4.9 


5.2 


LiH 


2.9 


8.3 


9.3 


6.8 


7.3 


H 2 


7.3 


19.1 


18.0 


14.8 


15.6 


CH 4 


11.6 


21.6 


20.2 


17.6 


18.2 


Na 2 


1.2 


4.6 


5.1 


4.4 


4.9 


Na 4 


0.4 


3.7 


4.2 


3.4 


3.9 



TABLE I: Calculated HOMO-LUMO gaps in electro-volts for 
several small molecules. Both the SCF and the NSCF ap- 
proximations produce gaps that are smaller and much closer 
to the HF and ASCF values than the LDA ones. The exper- 
imental values for LiH, Na 2 and Na 4 are 7.6, 4.5 and 3.4 eV, 
respectively. 
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FIG. 1: (color online) For the molecules shown, we compare 
the Hartree (H),Hartree-exchange (HX) contributions to the 
SCF gap correction (HXC) with the experimental (Exp), HF 
and LDA gaps (see text). All the lines are guides to the eye. 
The Hartree-only correction overestimates the experimental 
gap. Adding both exchange and correlation to the correc- 
tion lowers the gap value, bringing it closer to the HF and 
experimental values obtained from Ref. [34|- 



iment and are consistent with a simple screening picture. 
As argued above, these corrections are not included in the 
N-particle Kohn-Sham gap and hence are included in the 
so-called band-gap discontinuity, A xc , which is typically 
defined as the difference between the exact and Kohn- 
Sham gaps. In finite systems A^, should not be con- 
fused with the discontinuity in the XC potential, since, 
in addition to the discontuinity, there are continuous con- 
tributions originating for both the Hartree and XC po- 
tentials. In other words, the term band-gap discontinuity 
is not appropriate for finite systems. The discontinuity 
in the XC potential is not included in our calculations 
yet the effects we describe with Eq.© are not included 
in KS-DFT. 



C. Virtual Orbitals 

An interesting question is whether the application of 
Eqs.(0][n]) yields virtual orbitals in better agreement with 
quasiparticle orbitals. There is no precise way to evaluate 



4 



the quality of a mean field orbital, but since the adoption 
of Eqs.(g][5]) was based originally based on the property 
of the HF Hamiltonian discussed in the introduction, we 
take the HF LUMO as a reference state and compute 
its quantum-mechanical overlap with the LDA, SCF and 
NSCF LUMO's. In addition to this overlap we evaluate 
the following distance 

^£(l<f|-|cL,i|) 2 (6) 

where i labels a basis orbital, c^f is the i-th expansion 
coefficient of the HF LUMO, and is the corresponding 
coeffient of LDA, SCF and NSCF LUMO's. EqE] gives 
another measure of the similarity between the HF and 
LDA, SCF and NSCF LUMO expansion coefficients. The 
results of this comparison are shown in table [ill We note 
that the excellent agreement between the HF and the 
bare LDA LUMO's is improved in the molecules shown, 
except for the NSCF calculation of the sodium tetramer, 
where the HF LUMO overlaps with the NSCF LUMO+1 
instead of the the LUMO. 



Molecule 


Overlap 


Eq.© 


) 


LDA 


SCF 


NSCF 


LDA 


SCF 


NSCF 


H 2 


0.998 


1.000 


1.000 


0.050 


0.003 


0.012 


Li2 


0.975 


0.999 


1.000 


0.144 


0.050 


0.011 


LiH 


0.972 


0.997 


0.997 


0.122 


0.024 


0.026 


H 2 


0.994 


0.999 


1.000 


0.057 


0.019 


0.007 


CH 4 


0.997 


0.999 


0.996 


0.050 


0.016 


0.006 


Na 2 


0.955 


0.996 


0.999 


0.174 


0.056 


0.021 


Na 4 


0.948 


0.985 


0.000 


0.114 


0.058 


0.131 



TABLE II: Second, third and fourth columns: absolute val- 
ues of the overlaps between the Hartree-Fock LUMO and the 
LUMO's obtained by means of the LDA, SCF, and NSCF. 
Fifth, sixth, and seventh columns: distance (as given by 
Eq.©) between the HF and LDA, SCF, and NSCF LUMO 
states, respectively. 

While promising we do not expect these results to be 
general, due to the approximated nature of our calcu- 
lations. The XC functionals considered are continuous 
functionals of the density which are not self-interaction 
free. This has consequences both for the validity of our 
approach in infinite systems and for the values of the 
LUMO energy obtained for finite systems. Since the 



functionals we consider are continuous, one would ex- 
pect the correction represented by Eq.([5| to vanish for 
infinite systems. For the small molecules considered the 
electrostatic effect is large and its influence in the virtual 
orbitals seems to be captured by means of Eqs.([4][5|), un- 
fortunately this will not be the case for both extended 
systems or large molecules. The present approach does 
not modify the occupied orbitals, which would be iden- 
tical to those of underlying KS-DFT method, i.e., the 
calculated HOMO energies are too high because of the 
self-interaction error present in the functionals consid- 
ered. This leads to LUMO energies that are also too high 
when compared with experiments. In finite systems, the 
LUMO energy could in principle be improved by consid- 
ering functionals that are self- interaction free like, e.g., 
those of proposed in Refs. I35II361 . 



IV. CONCLUSIONS 

In conclusion, we have proposed a simple correction 
to the virtual orbitals of the KS-DFT Hamiltonian that 
approximates the screened electrostatic repulsion of the 
added particle by the other electrons. For the molecules 
tested, the method systematically improves on the calcu- 
lated fundamental gaps while provinding also improved 
LUMOs. Due to the various approximations used we do 
not regard the observed LUMO improvement as general. 
The non-self-consistent version of the proposed correc- 
tion improves over the KS-DFT LDA without a signifi- 
cant increase in computational time and difficulty of im- 
plementation. The proposed mean-field might be useful 
to improve the accuracy of KS-DFT-based calculations 
of the resistance of molecular junctions, since these have 
proven to be highly sensitive both to the gap and the 
details of the mean field orbital a ' 37 . Alternatively the 
proposed correction might be used to produce a better 
starting point for calculations of molecular systems based 
on many-body perturbation theory. 
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